#include "outline/grand_poly_bezier_builder.h"
#include "outline/bezier_guide.h"
#include <limits>

namespace outline {

template <typename T, typename R>
void grand_poly_bezier_builder<T, R>::reset(nodes_type const & polyline)
{
	reset_nodes(polyline);

	length_ = 3;
	estimations_.clear();
	estimations_.resize(size_);
}

template <typename T, typename R>
void grand_poly_bezier_builder<T, R>::set_params(grand_params<R> const & params)
{
	params_ = params;
}

template <typename T, typename R>
void grand_poly_bezier_builder<T, R>::inspect_greater_bezier_length()
{
	if (all_length_inspected())
		return;
	++length_;
	bezier_type b;
	for (uint i = 0; i < size_; ++i) {
		R e = bezier_guide(
				nodes_, i, i + length_,
				params_.low_max_iterations,
				params_.low_stop_mean_discrepancy,
				params_.low_oversample_factor,
				false, b
			);
		assert (e >= 0);
		estimations_[i].push_back(estimation_type(e, 
			b.start_control() - b.start(), b.finish_control() - b.finish()));
		assert(estimations_[i].size() == length_ - 3);
	}
}

template <typename T, typename R>
grand_poly_bezier_builder<T, R>::junction_layer::junction_layer(uint size, uint minl, uint maxl)
{
	assert (minl <= maxl);
	matrix.resize(size);
	l_min = minl;
	l_max = maxl;
	for (uint b = 0; b < size; ++b)
		matrix[b].resize(l_max - l_min + 1);
}

template <typename T, typename R>
R grand_poly_bezier_builder<T, R>::select_best(index_vector_type & ivec)
{
	if (!params_.low) {
		uint mj = std::min(params_.max_junctions, size_/3);
		ivec.resize(mj + 1);
		for (uint j = 0; j < mj; ++j) //uniformly distributed junctions
			ivec[j] = size_ * j / mj;
		ivec[mj] = ivec[0] + size_;
		return 0;
	}

	if (length_ <= 4)
		return -1;

	real_type best_D = std::numeric_limits<R>::max();
	uint best_b, best_j;
	if (length_ == size_ + 1) {
		for (uint b = 0; b < size_; ++b) {
			estimation_type const & e = inline_get_estimation(b, length_);
			real_type d = e.d;
			assert (d >= 0);
			if (params_.sharpness_penalty > 0)
				d += sharpness_penalty(params_.sharpness_penalty, e.d2, e.d1);
			if (d < best_D) {
				best_D = d;
				best_b = b;
				best_j = 0;
			}
		}
	}

	brick_type brick(params_.max_junctions);

	// copying the first layer
	junction_layer & jl0 = *(brick[0] = junction_layer_sptr(
		new junction_layer(
			size_,
			4, //min
			length_ //max
	)));
	for (uint b = 0; b < size_; ++b)
		for (uint l = jl0.l_min; l <= jl0.l_max; ++l)
			jl0(b, l) = step_type(l, estimations_[b][l-4].d);

	uint j;
	for (j = 1; j <= params_.max_junctions-1; ++j) {
		if (best_D <= params_.penalty_for_each_junction * j)
			break;
		junction_layer const & prev = *brick[j-1];
		junction_layer & curr =	
			*(brick[j] = junction_layer_sptr(new junction_layer(
				size_,
				prev.l_min + 3, //min
				std::min(prev.l_max + length_ - 1, size_ + 1) // max
			)));
		if (curr.l_min >= size_ + 1)
			break;
		for (uint b = 0; b < size_; ++b) {
			for (uint l = curr.l_min; l <= curr.l_max; ++l) {
				step_type & s = curr(b, l);
				s.null();
				uint m_min = std::max(4, (int)l - (int)prev.l_max + 1);
				assert (l > prev.l_min);
				uint m_max = std::min(l - prev.l_min + 1, length_);
				for (uint m = m_min; m <= m_max; ++m) {
					estimation_type const & e1 = inline_get_estimation(b, m);
					real_type d1 = e1.d;
					assert (d1 >= 0);
					uint b2 = (b+m-1) % size_;
					uint m2 = l-m+1;
					step_type const & second = prev(b2, m2);
					estimation_type const & e2 = inline_get_estimation(b2, second.m);
					real_type d2 = second.d;
					assert (d2 >= 0);
					real_type d = d1 + d2;
					if (params_.sharpness_penalty > 0) {
						d += sharpness_penalty(params_.sharpness_penalty, e1.d2, e2.d1);
						if (l == size_ + 1) { //fine the sharpness near the key node
							uint b0 = b2;
							uint l0 = m2;
							uint j0 = j-1;
							while (j0 > 0) {
								step_type const & s = (*brick[j0])(b0, l0);
								assert (!s.is_null());
								b0 = (b0 + s.m - 1) % size_;
								l0 -= s.m - 1;
								--j0;
							}
							assert ((b0 + l0 - 1) % size_ == b);
							estimation_type const & e0 = inline_get_estimation(b0, l0);
							d += sharpness_penalty(params_.sharpness_penalty, e0.d2, e1.d1);
						}
					}
					if (s.is_null() || s.d > d) {
						s.m = m;
						s.d = d;
					}
				}
				assert (!s.is_null());
				if (l == size_ + 1) {
					R D = s.d + j*params_.penalty_for_each_junction;
					if (D < best_D) {
						best_D = D;
						best_b = b;
						best_j = j;
					}
				}
			}
		}
	}

	ivec.clear();
	if (best_D == std::numeric_limits<R>::max())
		return -1;

	//back pass
	ivec.push_back(best_b);
	uint i = best_b;
	uint l = size_ + 1;
	j = best_j;
	while (j > 0) {
		step_type const & s = (*brick[j])(i % size_, l);
		assert (!s.is_null());
		i = i + s.m - 1;
		l -= s.m - 1;
		--j;
		ivec.push_back(i);
	}
	ivec.push_back(best_b + size_);

	return best_D;
}

template <typename T, typename R>
typename grand_poly_bezier_builder<T,R>::estimation_type const & 
grand_poly_bezier_builder<T,R>::get_estimation(uint begin, uint len, uint) const
{
	return inline_get_estimation(begin, len);
}

template <typename T, typename R>
void grand_poly_bezier_builder<T,R>::get_junction_types(
	index_vector_type const & ivec, std::vector<junc_type> & junc) const
{
	if (!params_.low) {
		for (uint i = 0; i != junc.size(); ++i)
			junc[i] = junc_cusp;
		return;
	}

	uint jcount = junc.size();
	assert(ivec.size() == jcount + 1);

	estimation_type const * e1 = &inline_get_estimation(ivec[jcount-1] % size_, ivec[jcount] - ivec[jcount-1] + 1);

	for (uint i = 0; i < jcount; ++i) {
		estimation_type const * e2 = &inline_get_estimation(ivec[i] % size_, ivec[i+1] - ivec[i] + 1);

		//if (square(e1->d2 + e2->d1) < 0.1 * (square(e1->d2) + square(e2->d1)))
		//	junc[i] = junc_symmetric;
		//else if (sharpness_penalty(real_type(1), e1->d2, e2->d1) < real_type(0.25)) //is angle between tangents less than 30 degrees
		//	junc[i] = junc_smooth;

		if (sharpness_penalty(real_type(1), e1->d2, e2->d1) < real_type(0.25)) //is angle between tangents less than 30 degrees
			junc[i] = junc_smooth;
		else
			junc[i] = junc_cusp;

		e1 = e2;
	}
}

template grand_poly_bezier_builder<unsigned short>;
template grand_poly_bezier_builder<short>;
template grand_poly_bezier_builder<int>;
template grand_poly_bezier_builder<float>;
template grand_poly_bezier_builder<double>;

} //namespace outline
